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ABSTRACT 

We investigate, by means of numerical simulations, the phenomenology of star for- 
mation triggered by low-velocity collisions between low-mass molecular clumps. The 
simulations are performed using an SPH code which satisfies the Jeans condition by 
invoking On-the-Fly Particle Splitting. 

Clumps are modelled as stable truncated (non-singular) isothermal, i.e. Bonnor- 
Ebert, spheres. Collisions are characterised by Mo (clump mass), b (offset parameter, 
i.e. ratio of impact parameter to clump radius), and M. (Mach Number, i.e. ratio 
of collision velocity to effective post-shock sound speed). The gas subscribes to a 
barotropic equation of state, which is intended to capture (i) the scaling of pre-collision 
internal velocity dispersion with clump mass, (ii) post-shock radiative cooling, and (iii) 
adiabatic heating in optically thick protostellar fragments. 

The efficiency of star formation is found to vary between 10% and 30% in the 
different collisions studied and it appears to increase with decreasing Mq, and/or 
decreasing b, and/or increasing M.. For b < 0.5 collisions produce shock compressed 
layers which fragment into filaments. Protostellar objects then condense out of the 
filaments and accrete from them. The resulting accretion rates are high, 1 to 5 x 
lO~ 5 M0yr -1 , for the first 1 to 3 x 10 4 yrs. The densities in the filaments, nn 2 > 
5 x 10 5 cm~ 3 , are sufficient that they could be mapped in NH3 or CS line radiation, 
in nearby star formation regions. 

Key words: hydrodynamics - method: numerical - stars: formation - fragmentation 
- binaries: general - ISM: clouds 



1 INTRODUCTION 

There is increasing observational evidence that cloud-cloud 
collisions account for a substantial fraction of the star for- 
mation in the Galaxy (Scoville, Sanders & Clemens 1986; 
Greaves & White 1991; Hasegawa et al. 1994; Miyawaki, 
Hayashi & Hasegawa 1999; Sato et al. 2000). Recent theoret- 
ical calculations suggest that such collisions provide a viable 
mechanism for triggering star formation (Tan 2000). Cloud- 
cloud collisions, like all dynamical star formation mecha- 
nisms, tend to result in the formation of groups of stars. 

Because the observational evidence for star formation 
triggered by cloud-cloud collisions comes principally from 
relatively violent star formation events that spawn massive 
OB stars, most numerical work to date has been concerned 
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with collisions between high-mass clouds. For example, the 
simulations of Chapman et al. (1992) treated collisions be- 
tween 75OM0 clouds. However, if substructure in molecular 
clouds is hierarchical with a low volume-filling factor (Scalo 
1985), one might expect a collision between two high-mass 
clouds to consist of many smaller-scale collisions between 
the lower-mass clumps of which the clouds are composed. 
Such a picture presupposes that the clouds and clumps are 
long-lived equilibrium entities. However, even in the cur- 
rently most popular paradigm of clouds being transient ob- 
jects formed and destroyed by turbulent motions (Padoan & 
Nordlund 2002; Mac Low & Klessen 2004; Bergin et al. 2004; 
Vazquez-Semadeni et al. 2005), fragmentation in shocks pro- 
duced by large-scale converging flows can be studied in terms 
of smaller-scale shocks triggered by collisions between equi- 
librium clumps provided that such structures are of suffi- 
ciently low-mass. 

The advantage of simulating collisions between lower- 
mass clumps is that the resolution requirements are less 
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severe. It has been shown, both in the context of Finite 
Difference simulations (Truelove et al. 1997; Truelove et al. 
1998), and in the context of Smoothed Particle Hydrody- 
namics (SPH) simulations (Bate & Burkert 1997; Kitsionas 
& Whitworth 2002), that gravitational fragmentation can 
only be modelled faithfully if the Jeans mass is resolved at 
all times. This requirement is normally referred to as the 
Jeans condition. Previous simulations of clump-clump col- 
lisions (Chapman et al. 1992; Pongracic et al. 1992; Turner 
et al. 1995; Whitworth et al. 1995; Bhattal et al. 1998; Mar- 
inho & Lepine 2000; Marinho, Andreazza & Lepine 2001) 
followed the evolution through several orders of magnitude 
in density and linear scale, and produced abundant high- 
mass protostellar fragments by gravitational fragmentation. 
However, these simulations did not satisfy the Jeans con- 
dition. Therefore, it is possible that real fragmentation on 
small scales was suppressed (Whitworth 1998), and/or that 
artificial fragmentation occurred. In this paper, we describe 
high-resolution simulations of clump-clump collisions, per- 
formed using SPH with On-the-Fly Particle Splitting (Kit- 
sionas & Whitworth 2002). Particle Splitting enables us to 
ensure that the Jeans condition is satisfied at all times. This 
is achieved with relatively modest computing resources, by 
introducing high resolution only in the regions where it is 
required. 

Two suites of simulations are presented aiming at i) the 
identification, free from numerical resolution constraints, of 
the dominant mechanism driving star formation in clump 
collisions and ii) the derivation of, at least, order of magni- 
tude estimates for the star formation efficiency in such col- 
lisions. In the first (and principal) suite of simulations, col- 
lisions between low-mass clumps having mass Mq = IOMq 
are simulated with various combinations of impact parame- 
ter and collision velocity. In the second suite, a few collisions 
between intermediate-mass clumps having Mo = 75M are 
simulated. The purpose of the second suite is twofold: (i) 
to explore the dependence on clump mass, and (ii) to re- 
peat some of the critical simulations reported by Bhattal et 
al. (1998) and establish which of the features they inferred 
might be attributable to inadequate resolution. 

In Section 2, we describe the physical model we use. 
In Section 3 we give a brief summary of our SPH code and 
explain how Particle Splitting is invoked to ensure the Jeans 
condition is always satisfied. A representative selection from 
the main suite of simulations (those using IOMq clumps) is 
presented in Section 4. Results from the second suite (75Mq 
clumps) are presented in Section 5 and compared with the 
results of Bhattal et al (1998). In Section 6 we discuss the 
results and summarise our main conclusions. 



2 PHYSICAL MODEL 

2.1 Clump and collision parameters 

In all cases, the collision involves two clumps of equal mass. 
In molecular clouds with hierarchical substructure, collisions 
between clumps at the same level of the hierarchy (i.e. with 
comparable mass) are the most probable (Scalo 1985) - al- 
though it would certainly be interesting to investigate colli- 
sions involving clumps of unequal mass. 

The collision is set up in the centre-of-mass frame, and 



so the clumps have antiparallel bulk velocities v c i ump and 
— v c iu m p. We define the Mach Number (M) of the collision 
as the ratio of the relative speed of the collision (2 |v c i ump |) 
to the effective post-shock sound speed (c s , see Eqn. 1 be- 
low). We define the offset parameter b as the ratio of the 
impact parameter of the collision to the clump radius. 

We take the clumps to be stable equilibrium isothermal 
spheres, i.e. they are contained by an external pressure and 
they are not singular (Ebert 1955; Bonnor 1956). In order 
to create such clumps we apply the procedure detailed in 
Appendix A4 of Turner et al. (1995). As demonstrated in 
Appendix B2 of Turner et al. (1995) this reproduces faith- 
fully the density profile of a stable isothermal sphere, with 
very little particle noise. 

In order to save computational time, the colliding 
clumps are touching when the simulations start. During the 
preceding approach, mutual tidal distortion should be small 
because the clumps move supersonically (M > 5). We define 
the a>axis to be parallel to the pre-collision velocities, and 
the y-axis to be parallel to the impact parameter. Hence the 
global angular momentum of the colliding clumps is parallel 
to the z-axis. 

In the context of a turbulent interstellar medium, the 
above model can be interpreted as the interaction of flows 
generated at the scale of the equilibrium clumps employed 
here either by sources driving turbulence at these scales or 
by the energy cascade of turbulence driven at larger scales 
(Padoan & Nordlund 2002; Mac Low & Klessen 2004; Bergin 
et al. 2004; Vazquez-Semadeni, Kim & Ballesteros-Paredes 
2005). The modelling of turbulent motions on smaller scales 
is not included self-consistently here, but rather it is repre- 
sented by an isotropic pressure force in the effective sound 
speed (Eqn. 1) of the clumps, in accordance with the scaling 
relations of Larson (1981) between the internal velocity dis- 
persion and the linear scale or the mass of a clump (e.g. see 
Eqn. 2). The use of an effective turbulent pressure to model 
turbulence acting on small scales suppresses the stochastic 
nature of any real density and velocity fluctuations expected 
from the presence of (supersonic or even transsonic) turbu- 
lent motions inside the clumps. Such an approach takes into 
account only the collective effect such fluctuations may have 
on the stability of the clump as a whole, neglecting the fact 
that small scale fluctuations can in turn induce collapse on 
these scales (Mac Low & Klessen 2004) , as we mainly want 
to identify the role of the collision itself on the triggering of 
star formation. 

Moreover, the current model does not take into account 
the effect that magnetic fields may have on the dynamics 
and the final outcome of the collisions. Although it is known 
that almost all studied clumps and cores are seen to be close 
to magnetic criticality (Crutcher et al. 2004), we effectively 
study here collisions between clumps threaded by weak mag- 
netic fields, which will not become dynamically important 
during the course of the collisions. Given the fact that i) 
collisions between low-mass clumps have not been studied 
before and ii) the numerical methods we employ provide 
an adequate framework for studying self-gravitating hydro- 
dynamics at the highest numerical resolution necessary, we 
believe that our assumption of magnetic subcriticality for 
the clumps is a reasonable first approach to this problem. 
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Figure 1. The equation of state given in Eqn. 1 presented 
as a density-temperature relation with cq given from Eqn. 2 
for the initial clump masses adopted here (see next session), 
c s ~ 0.2kms _1 , po = p c (pc is the initial central density of a 
clump; see next section), pi ~ 10~ 14 gcnr 3 . As soon as the col- 
lision starts the gas cools with T oc p~ 2 / 3 down to X ~ 10 K, a 
temperature that is reached at ~ 10~ 18 gcm~ 3 . Then it remains 
isothermal at T = 10 K, as it is thermally coupled to the dust, 
until it gets optically thick to its own cooling radiation and heats 
up adiabatically with T oc p 2 / 3 , for densities above 10 -14 gcm~ 3 . 



2.2 Equation of state 

We use a barotropic equation of state, 
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P ^ Po; 



P> Po, 



(1) 



that is illustrated in Fig. 1. The left hand side represents the 
square of the effective isothermal sound speed of the gas, i.e. 
when non-thermal pressure due to turbulence is included. 

At low densities, p ^ po (i.e. before the collision) the 
clump gas has an effective sound speed 
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in accordance with Larson's scaling relation between the 
mass and internal velocity dispersion of a clump (Larson 
1981). 

Once the density exceeds po = p c (p c is the initial cen- 
tral density of a clump; see the next section for the spe- 
cific values adopted here), i.e. essentially as soon as the gas 
is compressed by the collision shock, the gas cools down 
(with T oc p~ 2 / 3 , see the first term on the right hand side 
of Eqn. 1 for p > p ) to T ~ 10 K (Whitworth & Clarke 
1997; Whitworth, Boffin & Francis 1998). The post-shock 
sound speed is therefore c s ~ 0.2 kms -1 (corresponding to a 
cosmic mixture of H2 and He at 10 K). The gas reaches for 
the first time this temperature, and thus this sound speed, 



at ~ 10~ gem" and then it evolves isothermally for a few 
orders of magnitude in density (see Fig. 1). 

This part of the equation of state is similar to the 
piecewise polytropic equation of state advocated by Larson 
(2005). It differs mainly in the polytropic exponent assumed 
for the regime in which the gas is thermally coupled to the 
dust, i.e. above ~ 10~ 18 gem" 3 (Whitworth & Clarke 1997; 
Whitworth, Boffin & Francis 1998). For this regime, Larson's 
equation of state assumes slow heating of the gas (7 ~ 1.1), 
whereas we have taken here the gas to evolve isothermally 
(7 = 1). 

Finally, once the density in a collapsing protostellar ob- 



ject exceeds p± 



10- 



gcm , the gas is expected to be- 



come sufficiently optically thick to its own cooling radiation 
that it heats up adiabatically with T oc p 2 ^ 3 (Tohline 1982), 
given by the second term on the right hand side of Eqn. 1 
for p > po- We note that for a cosmic mixture of H2 and 
He at temperatures below 500 K, the adiabatic exponent is 
7 ~ 5/3 because the rotational degrees of freedom of H2 are 
frozen out; in our simulations the temperature does not rise 
above 200 K. 



2.3 The two suites of collisions 

A lOM© clump has equilibrium radius Ro — 0.22 pc, ef- 
fective isothermal sound speed Co — 0.35 km s~ (corre- 
sponding to a cosmic mixture of H2 and He at 35 K), 
central density p c ~ 2.6 x 10~ 20 gcm~ 3 and boundary 
density p h ~ 9.1 x 10~ 21 gcm~ 3 . With 1OM clumps 
we simulate collisions having all possible combinations 
of Mach Number M = 5, 10, and 15 (corresponding to 
|v c i ump | = 0.5, 1.0 and 1.5 km s _1 ) and offset parameter b — 
0.0, 0.2, 0.4, 0.6 and 0.8 (the parameters adopted for each of 
the simulations we performed are listed in the first seven 
columns of Table 1). The main results of this suite of simu- 
lations are presented in Section 4. 

A 75M clump has _Ro — 0.59 pc, Co ^ 0.60 km s -1 (cor- 
responding to a cosmic mixture of H2 and He at T ~ 100 K), 
p c ~ 1.1 x 10~ 21 gcm~ 3 and p h ~ 3.8 x 10~ 22 gcnT 3 . With 
75M© clumps we simulate only collisions having M = 9 
and 6 = 0.2, 0.4 and 0.5. We have chosen to study only this 
limited set of parameters for the intermediate-mass clump 
collision case, as a more detailed parameter study has been 
reported in Bhattal et al. (1998). Here we aim to verify with 
higher numerical resolution the validity of their conclusions 
on the phenomenology of clump collisions. In particular, we 
would like to investigate whether the existence of three dis- 
tinct mechanisms of star formation, which they identified in 
their clump collision simulations with varying b, depends on 
the low numerical resolution they employed. We have there- 
fore chosen to perform simulations with only three values 
of b, as according to Bhattal et al. (1998) collisions with 
b = 0.2, 0.5 and 0.4 produce protostars through rotational 
fragmentation, filament fragmentation and the combination 
of the previous two, respectively. In addition, we use this 
suite of simulations to investigate the effect of increasing 
clump mass on the outcome of clump collisions (by compar- 
ing the results of this suite of simulations with those of the 
low-mass clump collisions with M — 10 and corresponding 
values for b). The parameters adopted for each of the three 
intermediate-mass simulations we performed are also listed 
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Table 1. List of the initial condition parameters used for the complete set of simulations conducted in this paper (columns 1 to 7) as well 
as other characteristics of the evolution and final outcome of the simulations. Velocities are given in km s _1 , masses in Mq, temperatures 
in K, and times in Myr except for t cvo \ which is given in units of the tg ~ 0.05 Myr. 



in the first seven columns of Table 1. The main results of 
this suite of simulations are presented in Section 5. 



3 SPH WITH PARTICLE SPLITTING 

3.1 Standard self-gravitating SPH 

Smoothed Particle Hydrodynamics (SPH) (Gingold & Mon- 
aghan 1977; Lucy 1977; Monaghan 1992; Monaghan 2005), 
is a Lagrangian method for numerical hydrodynamics that 
assumes no symmetries or imposed grids, and is therefore 
very effective for treating problems that involve complex 3- 
dimensional geometries. SPH represents the fluid with an 
ensemble of TV discrete but extended particles. The parti- 
cles are overlapping, so that all intensive quantities can be 
treated as continuous functions both in time and space by 
averaging over neighbouring particles. To implement this, 
a smoothing kernel with compact support is used. The 
smoothing kernel describes the strength and extent of a par- 
ticle's influence. We use the polynomial M4 kernel (Mon- 
aghan & Lattanzio 1985). 

The SPH particles move with the fluid and all hydrody- 
namical properties are calculated at the particle positions. 
To evolve the ensemble of SPH particles, we use a system 
of three equations, namely the continuity equation, Euler's 
equation and a barotropic equation of state (Eqn. 1). The 
SPH formulations of the first two equations, giving the den- 
sity and acceleration of particle j, are 



Pi 



= E 



W 



(3) 



and 
dt 



E 



w 



+ 



' hi] 



Here r;j = Vj 



Vi and hij = 0.5(ft; + hj). W(s), W'(s) and 
W*(s) are the kernel, its derivative (W'(s) = dW/ds) and 
its volume integral (i.e. the mass fraction contained within 
s from the centre of the particle) . 

The first term on the right hand side of Eqn. 4 gives the 
contribution to the acceleration from hydrodynamic and ar- 
tificial viscosity forces. Artificial viscosity is included to pre- 
vent particle interpenetration at shocks. We use the stan- 
dard artificial viscosity (Monaghan 1992) 



where 



M»3 



o, 



Pij 



(vij Tij)<0; 
(Vij ■r ij )>0, 



|r l7 | 2 + 0.01ft?. ' 



(5) 



(6) 



v ij — v j ~ v i> Pij — 0.5(pi + pj) and Cij = 0.5(ci + cj) (the 
average isothermal sound speed). We have taken the artifi- 
cial viscosity parameters to be a = /3 = 1 . A test simulation 
of two flows colliding at Mach Number A^IO has shown 
that our code reproduces the expected density and velocity 
jump conditions of the shock. 

The second term on the right hand side of Eqn. 4 gives 
the gravitational contribution to the acceleration of particle 
j. By using the fact that the mass of a particle is kernel 
smoothed, and by invoking Gauss' gravitational theorem, 
close gravitational interactions (that would otherwise give 
unphysically large accelerations) are softened. It is implicit 
in Eqn. 4 that we use the same smoothing length to soften 
both gravitational and hydrodynamical forces. 

Due to the compact support of the kernel, the summa- 
tions in Eqn. 3 and the first term on the right hand side of 
Eqn. 4 are not over all particles, but over the small number 
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A/"ncib of nearby particles for which \vij \ < 2hij (= hi + hj). 
In three dimensions SPH gives good results with A/" ne ib ~ 50, 
and accordingly the smoothing length hi for particle i is ad- 
justed at each timestep so that jV nc ib ~ 50. 

Calculating the gravitational accelerations by a direct 
summation over all particle pairs, as implied by the right 
hand side of Eqn. 4, is an 0(Af 2 ) process, and therefore 
prohibitively expensive for large J\f. We avoid this expense 
by using the Tree-Code Gravity algorithm (Barnes & Hut 
1986; Hernquist & Katz 1989), which scales as O (TV log N). 
With this algorithm the computational domain is divided 
recursively into a hierarchy of cells within cells until, at the 
lowest level, each drawn cell contains either a single particle 
or no particle, from which we naturally store only tree-leaves 
containing a single particle each. Then we can calculate the 
gravitational interaction between particle j and particles at 
large distances from it, by treating these distant particles as 
a single point mass at the centre of mass of the cell to which 
they belong. 

For the time evolution of our code we use a second or- 
der Range-Kutta integration scheme and multiple particle 
timesteps. A detailed description of our standard code is 
given in Turner et al. (1995). 



3.2 The Jeans condition 

The Jeans condition requires that the Jeans mass Mjeans 
be resolved throughout the computational domain, at all 
times. In SPH, the minimum resolvable mass is estimated 
to be M rcso i ve( j ~ TVneibW-ptci, where m p t c i is the mass of 
a single SPH particle (Bate & Burkert 1997). In a three- 
dimensionally extended medium, the Jeans mass is given 
by Mjeans ~ G~ 3/2 p~ 1/2 c 3 . Thus the Jeans condition, 
A/resoived < Mjeans , reduces to a maximum density that can 
be resolved with SPH particles of mass m pt ci , 



that these fragments should have mass 



Mi 



frag 



(G 3 p clump A-!) 1/2 
/c s \ 3 M 
\c ) M 1 /' 2 



(10) 



where we have obtained the second expression by substi- 
tuting pdump ~ Cq/G 3 Mq for the pre-collision density of a 
clump in hydrostatic equilibrium, and the third expression 
by substituting for Co from Eqn. 2. If we now require that 
-^resolved < Mf rag , the maximum SPH particle mass becomes 



Afneib-M 1 / 2 



2 x 10~ 2 M ( 



o 



Mp 
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1/4 
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-1/2 



(11) 



Unless we consider very high- velocity collisions between very 
low-mass clumps, m^ax is much larger, and therefore less 
restrictive, than m max . 

It follows that, if the Jeans condition is to be satisfied, 
a standard SPH code cannot model a collision between two 
1OM clumps with fewer than ~ 400, 000 SPH particles, and 
a collision between two 75M© clumps requires ~ 3, 000, 000 
SPH particles. Simulations with such particle numbers are 
routinely performed at large parallel supercomputers (Bate, 
Bonncll & Bromm 2003; Bate & Bonnell 2005; Jappsen et al. 
2005), but on the smaller serial machines available to us 
they are prohibitive. In order to circumvent these prohibitive 
computational requirements, we have applied Particle Split- 
ting (Kitsionas & Whitworth 2002). 



G 3 (A/"neibm p tcl) 



(7) 



or equivalently a maximum SPH particle mass that can be 
used to model gas having density p, 



(8) 



G 3/2 p l/2_^ ne . b - 



With the barotropic equation of state that we are using 
(Eqn. 1), the combination c 3 / p 1 ^ 2 reaches its minimum when 
p ~ pi, i.e. when the gas switches from being approximately 
isothermal (c ~ c s ) for p < p < p\ to approximately adi- 
abatic for p > p\. Hence, the maximum SPH particle mass 
becomes 



G3/ 2 pl /2 A^neit 



5 x l0~ b M P 



© • 



(9) 



where the final evaluation has been made by substituting 
c s ~ 0.2 kms -1 , pi ~ 10~ 14 gem -3 , and A/" n eib ~ 50. 

A second - but as it turns out less severe - constraint 
on the mass of an SPH particle can be obtained by con- 
sidering the shock compressed layer which forms between 
two colliding clumps, and requiring that the fragments into 
which it breaks up are resolved. Whitworth et al. (1994a) 
have shown, on the basis of linear perturbation analysis, 



3.3 Particle Splitting 

In Particle Splitting, we replace individual particles (parent 
particles) with small groups of particles (families of child 
particles) either globally in a predefined sub-region of the 
computational domain where we anticipate the need for 
greater resolution (Nested Splitting) or conditionally accord- 
ing to some locally defined criterion (On-the-Fly Splitting). 
Each parent particle is replaced by thirteen child particles, 
having masses m c hiid = w- par ent/13 and smoothing lengths 
/'child = frijaront/13 1 / 3 . One of these child particles is placed 
at the same position as its parent, and the other twelve 
child particles are positioned on the vertices of an hexago- 
nal closed packed array, at equal distances £ = 1.5/i c hiid from 
the central child particle, and from each other. The family 
of thirteen child particles is then rotated to an arbitrary 
orientation. By positioning the child particles in this way, 
we ensure (i) that the family of child particles is as close as 
possible to the spherically symmetric mass distribution of 
the parent particle, and (ii) that the local transient fluctu- 
ations due to Particle Splitting are minimised (Kitsionas & 
Whitworth 2002). 

Each child particle is given a velocity v ch iid interpolated 
from the velocity field of its parent's neighbours i (including 
the parent j itself). This is formally given by 
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The initial densities and accelerations of the child parti- 
cles are then calculated using the standard SPH procedures 
(Eqns. 3 & 4). Subsequently, the child particles are also 
evolved with standard SPH procedures. The only difference 
is that, to mitigate interactions between adjacent particles 
having different masses, we have modified the scheme by 
which we calculate the smoothing lengths of particles. We 
now evolve hi, the smoothing length for particle i, so that 
the radius 2hi contains ~ 50 times the mass of particle i, 
i.e. 50m,i, rather than 50 other particles. This method can 
potentially introduce sampling errors in the calculation of 
the hydrodynamics quantities (see Eqn. 3 & 4), in case only 
a handful of neighbouring particles end up within 2hi, e.g. 
a child particle surrounded only by more massive particles. 
We have not addressed such special cases in our current code 
as corresponding tests, presented in Kitsionas (2000), have 
shown that such errors are negligible. However, we caution 
the reader that this can be a potential problem with the 
implementation of Particle Splitting in other SPH codes. 

In Kitsionas & Whitworth (2002) we have tested On- 
the-Fly Splitting on the standard Boss & Bodenheimer 
(1979) problem. On-the-Fly Particle Splitting is invoked in 
response to the imminent violation of the Jeans condition, 
i.e. whenever the density is about to exceed p m ax (Eqn. 
7). SPH with Particle Splitting (either Nested or On-the- 
Fly) reproduces faithfully the results obtained with a stan- 
dard high-resolution SPH simulation (using sufficient parti- 
cles to satisfy the Jeans condition without Particle Splitting, 
but also using much more memory and CPU), and the re- 
sults of high-resolution Adaptive-Mesh-Refinement Finite- 
Difference simulations (Truelove et al. 1997; Klein et al. 
1999). It also conforms to the analytic predictions of In- 
utsuka & Miyama (1992). 

In principle, Particle Splitting can be applied repeat- 
edly to produce successive generations of ever smaller SPH 
particles, and hence ever finer mass-resolution. However, in 
the simulations presented here only one generation of child 
particles is needed. By invoking a single generation of Parti- 
cle Splitting, we are able to start the simulations with SPH 
particles having m ptc i = 13m max ~ 6.5 x 1O _4 M0 (Eqn. 9). 
Therefore, initially a IOM0 clump can be modelled with just 
15,000 particles, and a 75M© clump with 110,000 particles. 
The critical density at which these initial particles have to 
be split to avoid violating the Jeans condition is 

Pspiit = 13~ 2 pi ~ 6 x 10~ 17 gcrrr 3 . (13) 

After this generation of Particle Splitting, adiabatic heating 
switches on before the simulation reaches its resolution limit 
again, and therefore, the Jeans condition is obeyed all the 
way up to the highest densities that can be achieved"!" 



t The Jeans mass increases with increasing density after adia- 
batic heating switches on, as Mj cans oc c 3 p _1 / 2 and in this regime 
c s oc p 1 / 3 (c/. Eqn. 1). We note that the Jeans mass increases with 
increasing density for all adiabatic exponents greater than 4/3. 



Figure 2. Linear density profile (particle plot) of one of the pro- 
tostars formed at the end of the simulation with Mo = IOMq, 
b = 0.2, Ai = 10, specifically the protostar at the bottom right 
of both panels of Fig. 4, with diameter of ~ 7.5 X 10~ 4 pc and 
density contrast of ~ 3 orders of magnitude. For comparison, the 
solid line illustrates the density profile of a single isolated par- 
ticle with mass and smoothing length corresponding to those of 
particles close to the peak density of the protostar. 

3.4 Plots and diagnostics 

We use grey-scale column density plots to present our re- 
sults. All structure formed is contained within a single layer 
and therefore such plots are not greatly confused by projec- 
tion effects. Column density plots are preferred to particle 
plots as the former give a more accurate representation of 
the total density field, and of what would be seen in opti- 
cally thin molecular-line (or dust-continuum) radiation, as- 
suming a uniform excitation temperature (or dust temper- 
ature). The figure captions give the linear size of each plot, 
the viewing axis, the time, and the range of the grey-scale 
(which is in all cases logarithmic, in units of g cm -2 ). 

In the sequel, we refer to all collapsing fragments iden- 
tified in our simulations as protostars and to their discs as 
protostellar discs. We would like to note that the terms "pro- 
tostar" and "protostellar" do not refer to pre-main-sequence 
(PMS) objects, as our simulations stop at much lower den- 
sity than that of PMS stars. Nevertheless, we use these terms 
for the fragments formed in our simulations, as all such frag- 
ments appear to have the characteristics of prestellar cores, 
i.e. collapsing cores not associated with IR sources. 

When a protostar forms, its linear density profile ap- 
pears like a normal distribution around its peak density. At 
the end of our simulations, this profile is rather steep (e.g. 
Fig. 2). To infer its mass, we take the protostar to extend 
out to di = 3o"i, where Oi is the FWHM in each of the three 
Cartesian axes i = x,y,z, and we sum the masses of the 
SPH particles that lie within all three diameters di. The 
radius of the protostar is given by the largest of the three 
radii r 4 = di/2. This method gives good results for the mass 
and the radius of spherical as well as disc-like protostars. 



High-resolution simulations of clump-clump collisions using SPH with Particle Splitting 7 



In essence, our method is similar to the Stutzki & Glisten 
(1990) clump-finding technique. Nevertheless, because of the 
low number of protostars formed in our simulations, we need 
not subtract each protostar from the density field before 
searching for the next one in the residuals. Therefore, we 
apply our method simultaneously to all protostars that are 

identified by eye in our simulations^ . Moreover, we note that, 
when applied to a few of our simulations, the Williams et 
al. (1994) clump-finding algorithm, as implemented for SPH 
by Klessen & Burkert (2000), gives results very similar to 
those of our method. In Table 1, we list the total mass of 
fragments, Mf rag , at the end of each simulation. 

Because our simulations evolve at very high resolution 
and due to the Lagrangian nature of SPH, the timestep be- 
comes extremely small in the densest parts of the computa- 
tional domain, i.e. within the protostellar fragments formed. 
As a result, after a few thousand timesteps our simula- 
tions approach a state of suspended animation, when the 
very small particle distances within the fragments dominate 
the system evolution. We have arbitrarily chosen to termi- 
nate all simulations after 10,000 timesteps as it becomes ex- 
tremely inefficient computationally to continue any further. 
To demonstrate that this is a reasonable choice, we have con- 
tinued one of the simulations for another 10,000 timesteps, 
during which the simulation advanced in physical time by 
only an additional ~ 1%. 



4 10M Q CLUMP COLLISIONS 

In these simulations the colliding clumps each have mass 
Mo — WMq, radius Ro — 0.22 pc, and pre-collision 
effective sound speed Co ~ 0.35 km s -1 . In the centre- 
of-mass frame, their velocities are v c i um p and — v c i ump 
and they collide with impact parameter bRo- We simu- 
late collisions with Mach Numbers M = 5, 10, 15 (cor- 
responding to |v c i ump | ~ 0.5, 1.0, 1.5kms _1 ), and b = 
0, 0.2, 0.4, 0.6, 0.8 (corresponding to impact parameters 
0, 0.044, 0.088, 0.132, 0.176 pc). 

Table 2 summarises the main results of all low-mass 
collision simulations performed. In the following subsections 
we describe in more detail a subset of the simulations. First 
we hold b = 0.2 (constant) and increase M through M = 

5 (§4.2), M = 10 (§4.3), and M = 15 (§4.4); then we hold 
M. = 5 (constant) and increase b through b — 0.2 (§4.2), b = 
0.4 (§4.5), and b = 0.6 (§4.6). In the following subsection 
we give a short overview of the expected outcome of clump 
collisions as discussed in Whitworth et al. (1994a; 1994b). 



4.1 Analysis 

The outcome of a collision between clumps can be analysed 
in terms of three factors, namely (i) the mass, extent and 
lifetime of the shock compressed layer, (ii) the fragmentation 
scale in the shock compressed layer, and (iii) the net angular 
momentum of the shock compressed layer. 



§ A protostar is identified for the first time when its peak den- 
sity exceeds ~100 times the density of its surroundings. We con- 
sider this time of first identification as the formation time of the 
protostar. 



(i) The mass, extent and lifetime of the shock com- 
pressed layer. For a given clump mass Mo, the amount of 
mass which is shock compressed depends on the impact pa- 
rameter. Evidently it is a maximum for head-on collisions 
(b — 0), and decreases monotonically with increasing b, be- 
coming negligible for b > 0.6. For larger impact parameter, 
two further factors come into play. Firstly, the mass of the 
shock compressed layer is increased somewhat at low Mach 
Number because the clump trajectories are focused by their 
mutual gravitational attraction. The escape speed from the 
surface of a clump, t> osc = (2GM a / Ro) 1 ^ 2 , is ~ 0.6 kms -1 
for a 10M Q clump, and so focusing is only important for 
the collisions with |v c i ump | ~ 0.5 kms -1 (i.e. M = 5). Sec- 
ondly, at high impact parameter and high Mach Number, 
a shock compressed layer forms but is then torn apart by 
shear before it can fragment. 

(ii) The fragmentation scale in the shock compressed 
layer. As discussed in Whitworth et al. (1994a; 1994b), the 
preferred fragment mass in a shock compressed layer is given 
by Eqn. 10, and the fragmentation scale (i.e. the mean sep- 
aration between fragments, in the plane of the shock com- 
pressed layer) is 



^frag 



Cs 



(Gp clurnp M) 1/2 

fcs\ Ro 
VcJ A4V2 

1/4 



M 

M 



Ro 



M 1 ' 2 



(14) 



We have obtained the second expression by substituting 
pciump ~ Co/G 3 Mo for the pre-collision density of a clump 
in hydrostatic equilibrium and Ro ~ GMo jc% for its radius, 
and the third expression by substituting for Co from Eqn. 
2. It follows that for given clump mass Mo and given (low) 
offset parameter b, increasing the Mach Number of the col- 
lision will decrease Lf rag , and hence increase the number of 
filaments into which the shock compressed layer fragments, 
and the number of protostars. At large b, the fragmenta- 
tion scale is increased somewhat by the shear in the shock 
compressed layer. However, this effect is only significant for 
b > 0.5, and the mass of the shock compressed layer is then 
negligible. 

We should stress that the fragmentation of a shock 
compressed layer is essentially a two-dimensional process, 
initially. The fragmentation scale Lf rag defines the size of 
proto-condensations in the plane of the layer. The thickness 
of the layer is much smaller than Lf rag . The number of fil- 
aments and protostars is therefore determined by the ratio 
of the lateral extent of the shock compressed layer to the 
fragmentation scale. 

(iii) The net angular momentum of the shock compressed 
layer. The net angular momentum of the shock compressed 
layer increases with increasing impact parameter and in- 
creasing collision velocity, i.e. increasing b and M. Increased 
angular momentum means that the ensemble of protostars 
has more collective orbital angular momentum, and so merg- 
ers between protostars are less likely. Additionally, it means 
that the material which accretes onto the individual proto- 
stellar discs from the filaments has higher specific angular 
momentum, and therefore spins them up more rapidly. 
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M\b 0.0 



0.2 



0.4 



0.6 



Two well-separated ro- 
tating objects. 0.68Mq 
in total. Only most mas- 
sive with spiral arms. 
Possible companions. No 
filaments. 



One spherical rotat- 
ing object. O.33M . No 
filaments. 



One disc-like object. 0.7 Mq. Spi- 
ral arms. No companions. A sin- 
gle filament. 



Two disc-like objects. 1.11M in 
total. Only most massive with 
spiral arms. Possible compan- 
ions. Single filament. 



10 One spherical rotating 

object. 0.35M Q . Two 
filaments. 



Two disc-like objects (+ a third 
forming). I.OMq in total. Both 
spiral arms. No companions. Net- 
work of filaments. 



Single disc-like object (+ a sec- 
ond forming). O.48M0. Spiral 
arms. Possible companions. Sin- 
gle filament. 



No shock. 



Single disc-like object. O.31M0. 
Spiral arms. No companions. No 
filaments. 



15 One disc-like rotating 

object. O.47M . No 
spiral arms. Network of 
filaments. 



Two disc-like objects (+ a 
third forming). O.95M0 in total. 
Only most massive with spiral 
arms. Possible companions. Well- 
defined network of filaments. 



No shock. 



Table 2. Summary of simulations and most important results for the low-mass clump collisions (Mo = IOJWq) for the different values 
of b and M. Shocks do not form in any of the b = 0.8 runs. 



In the following five subsections we describe five simula- 
tions which illustrate the above features of collisions between 
1OM clumps. 

4.2 M O = 1OM , 6 = 0.2, M = 5 

The initial conditions are shown in the left panel of Fig. 3. 
On-the-Fly Particle Splitting starts at t Bp nt ~ 0.463 Myr. A 
single tumbling filament forms at ~ 0.479 Myr, and mate- 
rial from the filament accretes onto a single central primary 
protostar T The filament is tumbling because the shock com- 
pressed layer from which it forms is tumbling, due to the or- 
bital angular momentum in a clump-clump collision at finite 
impact parameter. The steadily increasing specific angular 
momentum of the material accreting onto the protostar from 
the filament creates an accretion disc around the protostar, 
and spiral arms develop in this disc. We are unable to follow 
the simulation for long enough to determine whether these 
arms become sufficiently self-gravitating to condense out as 
secondary companions. 

At the end of the simulation (~ 0.496 Myr; right panel 
of Fig. 3), the mass of the protostar is ~ O.7OM and its ra- 
dius is ~ 90 AU. Its central density exceeds pi, which implies 
that it has started to heat up adiabatically. The number of 
active particles has increased from 30,000 to 55,300. 

4.3 M = 1OM , 6 = 0.2, M = 10 

The initial conditions are again as in the left panel of Fig. 
3. On-the-Fly Particle Splitting starts at t split ~ 0.339 Myr. 
By this stage, a network of tumbling filaments has started 
to form. At ~ 0.360 Myr, two protostars condense out of the 
filaments. The protostars are rapidly rotating, and accretion 
discs form around them shortly after their formation. The 
discs are approaching each other, and depending on their 

* We refer to protostars which form directly from the fragmen- 
tation of filaments as primary protostars. 



separation and their mutual alignment at periastron they are 
likely either to be captured into a binary, or to merge. There 
may be a third protostar starting to form a few timesteps 
before the end of the simulation. 

At the end of the simulation (~ 0.396 Myr; Fig. 4), the 
two protostellar discs have developed strong spiral arms, but 
they have not yet shown any inclination to fragment. The to- 
tal mass of the two protostars is ~ l.OM . The more massive 
protostar has mass O.59M and radius ~ 76 AU. The less 
massive protostar has mass O.41M and radius ~ 103 AU. 
Their central densities exceed pi, so they have started to 
heat up adiabatically. The minimum density in the filaments 
is pfii ~ 2.8 x 10~ 17 gem" 3 (n H2 ~ 5 x 10 6 cm" 3 ). The num- 
ber of active particles has increased from 30,000 to 64,300. 



4.4 M = 1OM , b = 0.2, M = 15 

Again, the initial conditions are as in the left panel of Fig. 
3. A network of tumbling filaments forms at ~ 0.306 Myr. 
When compared with the filaments formed in the previ- 
ous lower-velocity (Ai = 10) collision, the filaments of the 
M. = 15 collision are better defined, and they form on the 
whole surface of the shock compressed layer, not just at 
its centre. On-the-Fly Particle Splitting starts at t Bp \n ~ 
0.336 Myr. Two protostars condense out of the central fila- 
ment at ~ 0.348 Myr. These protostars are rapidly rotating, 
and accretion discs form around them shortly after their for- 
mation. The disc of the more massive (primary) protostar 
develops strong spiral arms. There are density enhancements 
at the points where the accretion flows intercept the spiral 
arms (right panel of Fig. 5). These density enhancements 
may subsequently become self-gravitating and condense out 
as secondary companions to the primary. There may be a 
third object starting to condense out in one of the other 
filaments. 

At the end of the simulation (~ 0.368 Myr; Fig. 5), the 
total mass of the two protostars is ~ O.95M . There is also 
another ~ O.1M associated with the spiral arms. The more 
massive protostar has mass O.53M and radius ~ 42 AU. 
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Figure 3. Column density plots for Mo = IOMq, b = 0.2, M = 5. Left Panel. Initial conditions viewed along the 
0.92 pc, Ay = 0.56 pc; the grey-scale is logarithmic, in units of g cm~ 2 , with sixteen equal intervals from 1.00 X 10~ a gcm~ z to 
2.69 X 10~ 2 gem -2 (or cquivalently, adopting solar composition, 2.50 X 10 2 ° H 2 cm~ 2 to 6.73 X 10 21 H 2 cm -2 ). Right Panel. The end of 
the simulation (t = 0.496 Myr) viewed along the z-axis; Ax = Ay = 0.016 pc; sixteen-interval logarithmic grey-scale, in units of g cm -2 , 
from 2.40 X lO" 1 gem" 2 to 2.04 X 10 3 gcm- 2 (6.00 X 10 22 H 2 cm -2 to 5.10 X 10 26 H 2 cm- 2 ). 



t=0.396 Myr 




-0.01 -5x10 3 5x10 3 

x (pc) 

Figure 4. Mo = IOMq, b = 0.2 M = 10 at the end of the simulation (t ~ 0.396 Myr). Left Panel. View along the z-axis; Aa; = 
0.016 pc, Ay = 0.014 pc; sixteen-interval logarithmic grey-scale, in units of g cm -2 , from 2.14 X 10 — 1 gem" 2 to 3.55 X 10 3 gcm — 2 
(5.35 X 10 22 H 2 cm- 2 to 8.88 X 10 26 H 2 cm- 2 ). Right Panel. View along the y-axis; Ax = Az = 0.02 pc; sixteen-interval logarithmic 
grey-scale, in units of gem -2 , from 1.05 X lO^gcm -2 to 1.78 X 10 3 gcm- 2 (2.63 X 10 22 H 2 cm~ 2 to 4.45 X 10 26 H 2 cm- 2 ). 
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Figure 5. Mo = 1OM , b = 0.2, M = 15 at the end of the simulation (t = 0.368 Myr), viewed along the rr-axis. Left Panel. Ay = Az = 
0.032 pc; sixteen-interval logarithmic grey-scale, in units of g cm , from 6.17 X 10 gem to 2.51 X 10 gem" - 2 (1.54 X 10 22 H 2 cm 
to 6.28 xlO 26 H 2 cm ). Right Panel. Zooming in on the protostar on the left edge of the left panel. Ay = Az = 0.008 pc; sixteen-interval 
logarithmic grey-scale, in units of g cm~ 2 , from 1.51 X 10" 1 gcm~ 2 to 1.62 x 10 3 gcm~ 2 (3.78 X 10 22 H 2 cm~ 2 to 4.05 x 10 26 H 2 cm" 2 ). 



The other protostar has mass O.42M and radius ~ 35 AU. 
Their central densities exceed pi, so they have started to 
heat up adiabatically. The minimum density in the filaments 
is pfu ~ 10~ 16 gcm~ 3 (nn 2 ~ 10 7 cm~ 3 ). The number of 
active particles has increased from 30,000 to 70,200. 



4.5 M = 1OM , 6 = 0.4, M = 5 

On-the-Fly Particle Splitting starts at t sp ii t ~ 0.525 Myr. A 
single tumbling filament forms at ~0.532 Myr. Two proto- 
stars form towards the two ends of the filament at ~0.535 
Myr. The protostars rotate rapidly, and accretion discs form 
around them shortly after their formation. One of the two 
protostellar discs develops strong spiral arms (right panel of 
Fig. 6), and there are density enhancements at the points 
where the spiral arms interact with the accretion flow. Sec- 
ondary companions to the primary protostar may subse- 
quently form from these enhancements. 

At the end of the simulation (~ 0.557 Myr; Fig. 6), the 
total mass of the two protostars is ~ 1.11M© (left panel 
of Fig. 6). The more massive protostar has mass O.67M0 
and radius ~ 77 AU. The other protostar has mass O.44M0 
and radius ~ 52 AU. Their central densities exceed pi , so 
they have started to heat up adiabatically. The protostars 
are separated by ~ 3500 AU. The minimum density in the 
filaments is p m ~ 2.43 x 10~ 17 gcm~ 3 (n H2 ~ 5 X 10 6 cm" 3 ). 
The number of active particles has increased from 30,000 to 
79,400. 



4.6 M = 10M©, 6 = 0.6, M = 5 

On-the-Fly Particle Splitting starts at i sp i it ~ 0.644 Myr. 
The two clumps move a long way into each other before 
the density increases significantly. Two single well-separated 
protostars form at ~ 0.678 Myr, but no filaments are formed 
in this collision. The protostars rotate rapidly, and accretion 
discs form around them shortly after their formation. One of 
the two protostellar discs develops strong spiral arms (bot- 
tom right protostar in Fig. 7). 

At the end of the simulation (~ 0.701 Myr; Fig. 7), the 
total mass of the two protostars is ~ O.68M0. The more 
massive protostar has mass O.4OM and radius ~ 53 AU. 
The other protostar has mass O.28M and radius ~ 85 AU. 
Their central densities exceed pi , and so they have started to 
heat up adiabatically. The separation between the two pro- 
tostars is > 30, 000 AU, and they are unbound. The number 
of active particles has increased from 30,000 to 50,000. 



4.7 Discussion of 10M© collisions 

As M is increased (with Mo = 1OM and 6 = 0.2 held con- 
stant), there are two main trends. Firstly, the fragmentation 
scale of the shock compressed layer, Lf ra g , becomes smaller, 
and therefore the layer breaks up into more filaments and 
forms more protostars. Secondly, the ensemble of protostars 
thus produced has more orbital angular momentum, and 
therefore mergers are less likely. 

As b is increased (with Mo = 1OM and M = 5 held 
constant), there are two main trends. Firstly, with larger 
impact parameter the mass of the shock compressed layer is 
reduced, and Lf rag is slightly increased. Consequently, fewer 
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Figure 6. Mo = IOMq, b = 0.4, Ai = 5 at the end of the simulation (t = 0.557 Myr), viewed along the z-axis. Left Panel. Ax = Ay = 
0.024 pc; sixteen-interval logarithmic grey-scale, in units of g cm - 2 , from 8.13 X lO^gcm" 2 to 2.14 X 10 3 gcm~ 2 (2.03 X 10 22 H 2 cm 2 to 
5.35 X 10 26 H2 cm -2 ). Right Panel. Zooming in on the protostar in the bottom right hand corner of the left panel. Ax = Ay = 0.01 pc; 
sixteen-interval logarithmic grey-scale, in units of g cm -2 , from 2.24 X 10 _1 gcm -2 to 1.58 X 10 3 gcm -2 (5.60 X 10 22 H2 cm -2 to 
3.95 x 10 26 H 2 cm -2 ). 




-0.05 0.05 
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Figure 7. M = IOMq, b = 0.6, M = 5, at the end of the 
simulation (t = 0.701 Myr), viewed along the z-axis; Ax = Ay = 
0.12 pc; sixteen-interval logarithmic grey-scale, in units of g cm -2 , 
from 2.45 X lO" 4 g cm" 2 to 2.75 X 10 2 g cm -2 (6. 13 X 10 19 H 2 cm" 2 
to 6.88 x 10 25 H 2 cm- 2 ). 



protostars are formed - and indeed for 6 > 0.6 no protostars 
are formed. Secondly, the pattern of fragmentation changes, 
in the sense that for higher impact parameter i) two pro- 
tostars form but with one being mainly material from one 
clump, and the other mainly material from the other clump; 
ii) the two protostars thus formed are not bound to one an- 
other. Fig. 8 summarises the above main conclusions. 

Accretion rates onto the protostars are typically 1 to 5 x 
Hn s M yr _1 . There is no obvious dependence of the accre- 
tion rate on the collision parameters. 

As we discussed in Section 2.2, the gas that gets com- 
pressed in the shock at the collision interface cools radia- 
tively down to ~ 10 K and then remains isothermal for a 
few orders of magnitude in density before it gets heated adi- 
abatically. Since gas fragmentation occurs in this isother- 
mal regime, we consider the free-fall time, ts, i.e. the most- 
commonly used measure of the fragmentation time-scale, to 
be that corresponding to gas of density 1 to 2 x 10~ ls g 
cm~ 3 , i.e the density at which gas reaches the isothermal 
temperature of 10 K for the first time (c/. Fig. 1). This 
gives tff ~ 0.05 Myr. Indeed, the time of fragmentation, tf lag 
(which is identical to the protostar formation time according 
to the definition given in Section 3.4), is in all simulations 
delayed by approximately one tg with respect to the time at 
which the gas in the simulation reached, for the first time, a 
density of ~ 10~ 18 g cm~ 3 . This means that there are two 
main phases of evolution in the simulations: in the first phase 
the dense shock compressed layer is forming (while the gas is 
cooling) and in the second (isothermal) phase the gas in this 
layer fragments. In Table 1, we list the fragmentation and 
the final time of each simulation, tf rag and t cn d respectively 
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Figure 8. Dependence of different quantities and phenomena 
on the increasing values of b and M for simulations of clump- 
clump collisions. Upward and downward pointing arrows indicate 
increasing and decreasing quantities, respectively. Note that the 
parameter space is divided in two sections: low-fe collisions pro- 
duce stronger shocks. Large-fe collisions reduce the clump inter- 
action. The transition happens at 6=0.4. 



(both given in Myr) , as well as the approximate duration of 
the second of the two phases defined above, which we call 
t cvo i (given in units of the iff), for which we have concluded 
that all collisions fragment exactly one tg ~ 0.05 Myr after 
the gas reached a temperature of 10 K for the first time. 

As can be seen from Table 1, the fragmentation time, 
ifrag, decreases with increasing Mach Number, Ai, and de- 
creasing offset parameter, b, as the formation of the shock 
compressed layer progresses faster with increasing clump ve- 
locity as well as in collisions where the clumps collide closer 
to head-on. 

In fact, the shock compressed layer breaks first into fil- 
aments which subsequently fragment to produce the (pri- 
mary) protostars. Fragmentation of the filaments is, there- 
fore, the main star formation mechanism at play in these 
low-mass clump collisions. It is the combination of gas cool- 
ing and self-gravity that are responsible for this behaviour. 
In particular, when we repeated the simulation with b — 0.2, 
Ai = 10 without taking into account the gas self-gravity, the 
two clumps quickly passed through each other as the shock 
compressed layer never became sufficiently dense (only by 
a factor of ~ 2.5 orders of magnitude) to break up into fil- 
aments and/or protostars. Repeating the same simulation, 
this time including the gas self-gravity but not allowing for 
gas cooling (i.e. the gas was isothermal at 35 K), no shock 
compressed layer formed but instead a single very massive 
(> 2.5M0) and highly bound object formed, which was ro- 
tating very fast and was attended by a disc with spiral arms. 

As can be seen from Tables 1 & 2, our simulations end 
at a stage when only about 1-5% of the total available gas 
mass has turned into protostars. Taking this value as a lower 
limit for the star formation efficiency (SFE) of the clump- 



clump collisions we model, we can also estimate an SFE 
upper limit from the gas mass that ends up gravitationally 
bound in the filaments. This gives a value of ~ 30-35% (e.g. 
for the 6 = 0.2, Ai — 10 run the bound mass in the filaments 
is ~ 6.7Mq). The early termination of the simulations does 
not allow us to measure the SFE in a more precise way than 
these estimated upper and lower limits. One can also use the 
fact that the protostars formed here show high mass accre- 
tion rates, characteristic of Class objects, and estimate an 
expected final mass for the protostars ~ 1 Myr after their 
formation (using mass accretion rates lower by one or two 

orders of magnitude for the duration of the Class I phased). 
In this way one ends up with a final mass of 1 — 2Mq per pri- 
mary protostar formed. Using also the fact that 1-2 primary 
protostars form on average in our simulations**, we infer 
a SFE of order 10-20% for the low-mass collisions. Based 
on the number of Jeans masses in the filaments, we expect 
that in the high Ai, b — 0.2 collisions a third primary pro- 
tostar will also form, as the number of Jeans masses in the 
filaments at the end of these collisions is of order 3. Thus, 
we expect the SFE of such collisions to be higher than the 
average, but still in the 10-20% range. 

In some of our simulations, spiral arms form in the pro- 
tostellar discs and we find some evidence for interactions 
between the spiral arms and the accretion flows. However, 
our simulations stop at an early stage of the disc evolu- 
tion due to timestep constraints. Thus, we can not confirm 
that such interactions are efficient in forming companions 
to the protostars. Formation of secondaries by accretion- 
induced rotational instabilities and/or disc-disc interactions 
would increase the SFE. We note, however, that the Toomre 
q parameter^ for most of the discs (e.g. both discs of the 
6 = 0.2, Ai — 10 collision; see Fig. 4) is steadily decreasing 
and it is q > 1 at the end of the simulations. 



4.8 Repeating the M = 1OM , b = 0.2, Ai = 10 
simulation with sink particles 

4-8.1 The algorithm and sink parameters 

In order to derive a more direct estimate of the number 
and the final mass of the protostars that form in such colli- 
sions, we have repeated the b = 0.2, Ai = 10 collision this 
time including sink particles (Bate, Bonnell & Price 1995). 
In particular, if certain conditions are satisfied (protostellar 
peak density exceeding a certain density threshold while the 
protostar is gravitationally bound; for details see below), we 
force all gas particles within a protostar to be removed from 
the simulation and replaced with a single collisionless (star) 
sink particle having as mass the sum of the masses of the 



II We assume, very crudely, that the Class phase lasts for ~ 0.1 
Myr and the Class I phase for ~ 1 Myr, after which the mass 
accretion rate becomes negligible. 

** The number of Jeans masses in the filaments is on average of 
order 2 at the end of the simulations. 

ft To calculate q, we use the average values of the sound speed, 
density, velocity, distance from the disc centre, and distance from 
the disc midplane for all particles in a disc. 
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Figure 9. Continuing the Mo = WMq, b = 0.2 M. = 10 simulation by using sink particles (overlayed solid circles; note that the 
symbol used for the sinks is larger than the actual sink radius at the scale of these plots). View along the y-axis; Ax = Az = 
0.02 pc in all panels. Upper Left Panel, t ~ 0.393 Myr, i.e. comparable to the right panel of Fig. 4 (the end of the corresponding 
cirvinict^n •'"'fii««( ointoV ;vt mn _mtor,Bii logarithmic grey-scale, in units of g cm -2 , from 3.31 X 10 — 1 gcm -2 to 5.25 X 10 2 gcm~ 2 

i~ 2 ). Upper Right Panel, t ~ 0.402 Myr; sixteen-interval logarithmic grey-scale, in units of g 
cm~~, trom 'A.titi X 1U~" gem - " to 2.04 X 10 gcm~ 2 (7.20 X 10 22 H 2 cm~ 2 to 5.10 X 10 25 H 2 cm~ 2 ). Both Lower Panels, t ~ 0.417Myr; 
the Bottom Right Panel is an exact copy of the Bottom Left Panel after having shifted the latter by 0.008 pc along the z-axis (in order 
to allow for direct comparisons with the corresponding plots of Fig. 10); sixteen-interval logarithmic grey-scale, in units of g cm -2 , from 
2.45 X 10" 1 gem" 2 to 9.12 X 10 1 gem" 2 (6.13 X 10 22 H 2 cm~ 2 to 2.28 X 10 25 H 2 cm' 2 ). 
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1.28 x 10 22 H 2 cm 
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to 1.31 X 10 26 H 2 cm- 2 ). 



gas particles it replaces"". The position and the velocity of 
a sink are given by mass- weighted sums of the positions and 



M It is essential to introduce sinks automatically in the simula- 
tions, as the replacement by hand of protostars with sinks will 
lead the system to a new state of suspended animation as soon 



the velocities, respectively, of the gas particles that the sink 
replaces. Subsequent to their formation, sinks are allowed 
to accrete more gas that enters their radius of influence, 



as another protostar forms, in which case we will have to replace 
this new protostar with another sink and so on. 
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Figure 10. Continuing the Mo = 10M©, b = 0.2 M = 10 simulation by using sink particles (overlayed solid circles; note that the 
symbol used for the sinks is larger than the actual sink radius at the scale of these plots). View along the j/-axis; Ax = Az = 0.02 pc 
in both panels (to be compared directly with the bottom right panel of Fig. 9). Left Panel, t ~ 0.435 Myr; sixtccn-interval logarithmic 
grey-scale, in units of g cm~ 2 , from 3.09 X 10 _1 gcm~ 2 to 9.55 X 10 2 gcm~ 2 (7.73 X 10 22 H 2 cm -2 to 2.39 X 10 26 H 2 cm -2 ). Right Panel. 



At the end of the simulation (t ~ 0.440 Myr); sixtecn-intcrval logarithmic grey-scale, in units of g cm 
6.61 x 10 2 gcm- 2 (8.48 x 10 22 H 2 cm~ 2 to 1.65 x 10 26 H 2 c m - 2 ). 
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Figure 11. Left Panel. The positions and velocities of the 9 sinks shown in the right panel of Fig. 10, i.e. at the end of the simulation 
with sinks. Circles denote the binary components (M > IMq), squares denote low-mass protostars (O.IA/q < M < IMq) and triangles 
denote protostars in the sub-stellar mass range (M < O.IMq). The lines indicate the projection of the velocity of each sink on the 
xz-plane using a scaling of 1000 km s~ x : 1 pc. A 1 km s^ 1 velocity line is indicated for clarity at the bottom of this panel. The direction 
of each velocity vector is pointing away from the corresponding symbols. Right Panel. Zooming in close to the binary (the region enclosed 
by the square in the left panel). For clarity the velocity of the low-mass companion of the binary is illustrated by a dashed line. The 1 
km s^ 1 velocity line is indicated at the top of this panel. 
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r S ink- This is the only way the sinks interact with the gas, 
apart from the boundary conditions that they impose on 
their neighbouring gas particles (see the discussion below). 
Gas particles that are accreted by a sink are removed from 
the simulation. The sinks are subject only to the force of 
gravity, exchanging gravitational interactions with all types 
of particles. Sink-sink interactions are Plummer softened at 
distances smaller than 2 x r s i n k- 

The conditions for sink creation are that a protostar is i) 
dense enough, i.e. its peak density exceeds a certain thresh- 
old, specifically p s ink = 10~ 12 g cm~ 3 , which is 100 times 
higher than pi, the density after which adiabatic heating 
switches on; and ii) gravitationally bound, i.e. the gas within 
the protostar has already started collapsing. In terms of the 
Jeans analysis, from the second of the above sink formation 
conditions we obtain the minimum mass of a sink that can 
be resolved at p s i n k as well as the corresponding sink cre- 
ation radius which also serves as the sink accretion radius, 
r s i n k. In particular, using the fact that for densities between 
pi = 10~ 14 and p s ink = 10~ 12 g cm~ 3 the gas evolves as 
c oc p 1 / 3 (cf. Eqn. 1) the minimum fragment mass at pi 
of Eqn. 9 translates to a minimum mass of a sink at p s i n k 
of order 0.05M Q for A/" nc ib ~ 100 (Bate & Burkert 1997). 
Therefore, we obtain r s i n k ~ 20 AU. 

Gas particles are accreted by a sink when i) they have 
approached the sink at a distance less than r s i n k, and ii) they 
are gravitationally bound to the sink. In order to model ac- 
cretion by sinks properly, a number of boundary corrections 
must be applied to the neighbouring gas particles of a sink, 
so that these gas particles do not feel a discontinuity in the 
density field due to the proximity of a sink. We have im- 
plemented all the boundary corrections given by Bate el al. 
(1995). 

4.8.2 Results 

The two protostars identified at the end of the simulation 
without sinks (Fig. 4), are replaced by sinks at t ~ 0.391 and 
t ~ 0.393 Myr (upper left panel of Fig. 9; these are the sinks 
on the right and in the middle of the panel, respectively). 
Note that the sink on the right hand side of the upper left 
panel of Fig. 9 lies below (lower y value) the sink in the 
middle of the panel (cf. the left panel of Fig. 4). As the 
simulation with sinks progresses, the sink on the right hand 
side of the upper left panel of Fig. 9 first moves to the left 
towards the centre of the computational domain and later 
it is forced to move along the z-axis by shear produced by 
material entering the shock compressed layer (at an angle) 
from the opposite side. At the time of the upper right panel 
of Fig. 9, this sink has exited the small field of view of this 
panel (from the top side). 

In the meantime, a third protostar formed (cf. Section 
4.3) in the top left hand corner of the displayed region and 
at t ~ 0.400 Myr it was replaced by a sink. At the time of the 
upper right panel of Fig. 9 (t ~ 0.402 Myr), the two sinks 
remaining in the displayed region have already approached 
each other along the filament. They have a close encounter 
with periastron ~ 200 AU at t ~ 0.403 Myr and capture 
themselves into a binary. A large circumbinary disc forms 
around the binary. The binary components continue accret- 
ing matter from this disc and their orbit is hardened. At 
the same time, the disc material is replenished by material 



falling on it from the filaments. Interactions between the ac- 
cretion flows and the disc lead to the formation of further 
companions to the binary (see e.g. the sinks on the spiral 
arms of the disc in the lower panels of Fig. 9). Most of these 
companions are ejected within a few orbits around the bi- 
nary after having a close encounter with it (e.g the sinks in 
the lower density regions of the left panel of Fig. 10). 

At the end of the simulation with sinks (right panel of 
Fig. 10, t ovo i ~ 2.6 iff), the binary components have mass 
of ~ 1.5Mq each, and their separation is ~ 40 AU, i.e. the 
minimum separation resolved due to the Plummer softening 
that is used to model sink-sink interactions. Eight additional 
sinks have formed in total (including a sink that has exited 
the displayed region; this sink has final mass of ~ 0.5M©). 
The masses of the secondary protostars, grouped in two mass 
bins, are denoted by different symbols in Fig. 11 (squares for 
0.1M© < M < 1.0M© and triangles for M < 0.1M©). Their 
projected velocities are also indicated by the lines in Fig. 
11. None of them remains bound to the binary. However, 
they have been ejected from the system with moderate ve- 
locities of order a few km s _1 (e.g. the object having a close 
encounter with the binary on the right panel of Fig. 11). 

A more detailed account of the properties of this sys- 
tem of protostars is beyond the scope of this paper. We 
note, however, that the general picture drawn from our sim- 
ulation with sinks with respect to the formation and evolu- 
tion of this system of protostars is similar to that of one of 
the subclusters reported by Bate & Bonnell (2005), which 
has also formed along filaments. A small difference to the 
results of Bate & Bonnell is the lower rate with which sec- 
ondaries form in our simulation, which we attribute to the 
stiffcr equation of state used here (with polytropic exponent 
of 5/3 instead of 7/5 that is used by Bate & Bonnell). A 
stiff equation of state provides proto-fragments with larger 
thermal pressure support. For example, we have noted that 
fragments like the one in the top left hand corner of the 
right panel of Fig. 10 form a sink only when they fall onto 
the disc, i.e. when they get very quickly loaded with mass 
exceeding a Jeans mass. We also note that no individual 
discs and/or envelopes can be resolved here for any of the 
10 sinks formed, i.e. if any of these protostars is attended by 
a disc and/or surrounded by an envelope, this disc/envelope 
will be smaller than r sink ~ 20 AU in size. The separation 
of the binary is, accordingly, expected to be smaller than 
~ 40 AU which is the minimum separation than can be re- 
solved. Finally, we note that, despite the very small number 
statistics, the mass distribution of the sinks formed here is 
consistent with the observed protostellar mass distributions 
in several star forming regions. 

The total mass in sinks at the end of the simulation is 
5.1 Mq, which implies a SFE of order 25%. There is, how- 
ever, about 3 Mq more bound gas in the filaments, but we 
expect only the binary components to have access to this 
material. Moreover, we expect accretion onto these two pro- 
tostars to start decreasing at subsequent times, i.e. when 
they enter the Class I phase, during which stellar feedback 
effects (which we do not model here) start affecting the disc 
from which they now accrete. We do not anticipate, there- 
fore, more than 1-2 Mq of additional stellar mass in the bi- 
nary and/or additional secondaries, limiting the SFE of this 
collision to a maximum value of ~ 30%. Because this colli- 
sion is one of the most efficient among the collisions we study 
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here, we estimate the average SFE of the low-mass clump 
collisions to be of order 20-30%. After having conducted the 
simulation with sinks, we note that, in terms of primary pro- 
tostellar mass, our initial SFE estimate, which was made in 
the previous section by extrapolating the masses obtained 
from the simulations without sinks and using well known 
protostellar mass accretion rates and lifetimes, was rather 
accurate. The SFE estimate obtained from the simulation 
with sinks also accounts for the formation of secondary pro- 
tostars, which we could not follow before introducing sinks, 
and has lead to this higher SFE of order 20-30%. 



5 75M Q CLUMP COLLISIONS 

In this section we present the results of the suite of three 
simulations involving collisions between clumps having mass 
M = 75M Q . The Mach Number is set to M = 9, and 
the offset parameter is set to b — 0.2, 0.4, and 0.5. The 
purpose of this suite of simulations is to explore how the 
results depend on clump mass Mo , and to test the reliability 
of the results reported by Bhattal et al. (1998). 

The main effect of increasing the clump mass, Mo, and 
hence also the clump radius, Ro, is that - at fixed b - the 
shock compressed layer is more extended and the fragmen- 
tation scale is also somewhat greater (c/. Eqn. 14). Conse- 
quently, the shock compressed layer breaks up into a more 
complex network of filaments, as the number of filaments 
should be tx M^ 4 . 

The main conclusion concerning the simulations re- 
ported by Bhattal et al. (1998) is that their results are seri- 
ously corrupted, both by the fact that the Jeans condition 
is violated, and by the fixed (and rather large) gravity soft- 
ening length. 
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Figure 13. M = 75M e , b = 0.2, M = 9, at the end of 

the simulation, (t = 0.640 Myr), viewed along the y-axis to 
show the network of filaments produced; Ax = Az = 0.2 pc; 
sixteen-interval logarithmic grey-scale, in units of g cm -2 , from 
7.24 X 10- 2 gem" 2 to 1.51 X 10 2 gem" 2 (1.81 X 10 22 H 2 cm" 2 to 
3.78 x 10 25 H 2 cm- 2 ). 



5.1 M = 75M , 6 = 0.2, M = 9 

After ~ 0.61 Myr, the shock compressed layer has become 
sufficiently dense and massive to break up into a network of 
tumbling filaments. Protostars then start to condense out at 
the intercepts of the filaments. On-the-Fly splitting starts at 
~ 0.63 Myr. In total, four protostars form, but two of them 
quickly merge. The three surviving protostars are rotating 
rapidly, and attended by accretion discs with strong spiral 
structure. Because of the increasing offset between the op- 
posing streams accreting onto a protostar along the filament, 
the protostellar disc is spun up, and we anticipate that the 
disc will eventually fragment to produce secondary compan- 
ions. However, the simulation has been terminated before 
this happens. The accretion streams are very smooth. 

In comparison with the lower-mass (IOMq) collision at 
b = 0.2, more filaments and more protostars are formed, 
because the area of the shock compressed layer is larger and 
the fragmentation scale Lf rag is only moderately larger. 

At the end of the simulation (~ 0.64 Myr; Figs. 12 & 
13), the total mass of the three protostars is ~ O.85M0, and 
their radii are in the range 100 to 130 AU. Their central den- 
sities exceed pi, so they have already started to heat up adia- 
batically. The separation between the protostars is ~ 10 4 AU 
(Fig. 12), and they are only weakly bound, so it is almost 
certain that within 1 Myr one of them will be ejected grav- 
itationally by the other two. The minimum density in the 



filaments is pm ~ 3.0 x 10 17 gem 3 (tih 2 ~5x 10 6 cm 3 , 
Fig. 13). The number of active particles has increased from 
220,000 to 260,000. 

In contrast, Bhattal et al. (1998) found that for low-6 
collisions a single primary protostar formed at the centre of 
the collision and accreted material from a single filament. 
The primary protostar rapidly acquired an accretion disc, 
and grew in mass and angular momentum as the offset be- 
tween the opposing accretion streams along the filament in- 
creased. Due to the development of spiral structure in the 
disc, and the lumpiness of the accretion streams, lower-mass 
secondary companions formed in the disc around the pri- 
mary. There was no suggestion of the shock compressed 
layer fragmenting into multiple filaments. Moreover, the sec- 
ondary companions modelled by Bhattal et al. were not al- 
ways resolved properly, as in some cases they contained fewer 
than 50 particles. Also Bhattal et al. used a large and con- 
stant gravity softening length e ~500 AU, so that they could 
not resolve the formation of an object until its size became 
>500 AU. 

We conclude that the evolution of the density field is 
followed more faithfully in our current simulations with Par- 
ticle Splitting. The fact that the Jeans condition is obeyed, 
prevents a central object from forming artificially before the 
filaments. Also there is no preferred length scale, so we can 
trust the detailed evolution of the discs and their dynamical 
interaction with the filaments. 
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Figure 12. M = 75M Q , b = 0.2, M = 9, at the end of the simulation, (t = 0.640 Myr). Left Panel. View along the z-axis; Ax = 0.04 pc, 
Ay = 0.08 pc; sixteen-interval logarithmic grey-scale, in units of g cm~ 2 , from 3.63 X 10~ 2 gcm~ 2 to 7.59 xlO 2 gcm~ 2 (9.08xl0 21 H2 cm~ 2 
to 1.90 X 10 26 H2 cm -2 ). Right Panel. View along the y-axis; Ax = Az = 0.028 pc; sixteen-interval logarithmic grey-scale, in units of g 
cm -2 , from 2.00 X 10" 1 gem" 2 to 1.70 X 10 3 gem" 2 (5.00 X 10 22 H 2 cm" 2 to 4.25 X 10 26 H 2 cm- 2 ). 



5.2 Mo = 75M Q , 6 = 0.4, M = 9 

A single tumbling filament forms at ~ 0.66 Myr. On-the-Fly 
Particle Splitting starts at about the same time. At least 
four protostars condense out of the filament; a further two 
may have started to condense out. 

At the end of the simulation (~ 0.686 Myr; Fig. 14), 
the total mass of the four well established protostars is ~ 
O.4OM0. They are positioned randomly along the filament, 
and they are falling towards one another. The most massive 
protostar is disc- like with a radius of ~ 170 AU (Fig. 14). 
The three smaller protostars are still roughly spherical with 
radii ~ 45 AU. Their central densities exceed pi, so they 
have started heating up adiabatically. The number of active 
particles has increased from 220,000 to 260,000. 

This is very similar to what happened in the corre- 
sponding simulation performed by Bhattal et al. (1998), ex- 
cept that their simulation was followed further, and the disc 
round the most massive protostar was therefore spun up 
by accretion along the tumbling filament. As a consequence 
it became rotationally unstable, and fragmented to produce 
secondary companions. This might also have happened here, 
if we had been able to follow the simulation for long enough. 
However, once again we can have little faith in the Bhattal et 
al. result because some of the secondaries contained at their 
inception less than 50 particles, and gravity was severely 
softened on scales < 500 AU. 

5.3 M = 75M Q , 6 = 0.5, M = 9 

On-the-Fly Particle Splitting starts at ~ 0.72 Myr. At ~ 
0.73 Myr, a single tumbling filament forms, and fragments 



into several protostars. Two of them, at the bottom right 
hand corner of Fig. 15, condense out first and eventually, 
at ~ 0.746 Myr, they merge. Subsequently, another three 
protostars form, two near the centre of the filament and one 
at the left hand end (Fig. 15). All fragments remain spherical 
and are strongly centrally condensed. Some other density 
enhancements appear in the filament, suggesting that more 
fragments may be forming. 

At the end of the simulation (~ 0.75 Myr), the total 
mass of the four well established protostars is ~ 0.65Mq. 
The most massive one has radius ~ 35 AU and the radii of 
the other three are ~ 100 AU. Their central densities ex- 
ceed pi, so they have started to heat up adiabatically. The 
minimum density in the filament is pm ~ 2 x 10~ 18 gcm~ 3 
(iih 2 ~5x 10 cm ). The number of active particles has 
increased from 220,000 to 342,000. 

Again, this is very similar to what happened in the cor- 
responding simulation performed by Bhattal et al. (1998), 
where a single tumbling filament formed and then frag- 
mented. Bhattal et al. performed a standard resolution sim- 
ulation of this case, which produced two protostars after 
~ 1 Myr, and a higher-resolution simulation in which each 
component of the binary was itself a binary, i.e. an hierar- 
chical quadruple. 

5.4 Discussion of 75M Q collisions 

In the b = 0.2 collision, the shock compressed layer has 
considerable lateral extent, and fragments into a network of 
tumbling filaments; multiple protostars then condense out of 
these filaments. In the b = 0.4 and b = 0.5 collisions, the ex- 
tent of the shock compressed layer is considerably smaller, 
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Figure 14. M = 75M Q , b = 0.4, M = 9, at the end 
of the simulation (t = 0.686 Myr), viewed along the z-axis; 
Ax = Ay = 0.008 pc; sixteen-interval logarithmic grey-scale, in 
units of g cm- 2 , from 9.33 X lO" 1 gem" 2 to 7.41 X 10 2 gem" 2 
(2.33 x 10 23 H 2 cm- 2 to 1.85 x 10 26 H 2 cm' 2 ). 



t=0.746 Myr 




Figure 15. M = 75M , b = 0.5, M = 9 at the end 
of the simulation (t = 0.750 Myr), viewed along the z-axis; 
Ax = Ay = 0.008 pc; sixteen-intervals logarithmic grey-scale, 
in units of g cm -2 , from 3.98 X 10 -1 gem -2 to 6.03 X 10 3 gem -2 
(9.95 x 10 22 H 2 cm" 2 to 1.51 x 10 27 H 2 cm' 2 ). 



and in the first instance it produces a single tumbling fil- 
ament, which then fragments into a line of protostars. As 
these protostars accrete material with increasing specific an- 
gular momentum from the tumbling filament, accretion discs 
form around them and develop spiral structure. At the same 
time the protostars fall towards one another, and they are 
therefore likely to interact, either merging, being captured 
into binary systems, or triggering the formation of further 
protostars. 

The main difference between the IOMq collisions of Sec- 
tion 4, and the 75Mq collisions of this section, is that the 
latter involve more mass and produce a shock compressed 
layer of greater extent. Consequently the filaments are longer 
and spawn more protostars, typically 4 to 6 (as compared 
with 1 to 2 for the IOMq collisions). Despite the fragmen- 
tation scale being somewhat larger, there tend to be more 
filaments, at least for low-& collisions. 

Taking into account that in the 75Mq collisions i) there 
are on average three times more primary protostars than in 
the IOMq collisions, and ii) there is more mass available 
for protostars to accrete in the shock compressed layer (and 
thus also in the filaments), and using the SFE estimated in 
Section 4.8 for the IOMq collisions, we arrive to an average 
SFE estimate of order 10-15% for the 75Mq collisions (as 
the total available gas mass is 7.5 times larger). 

It seems inescapable that the results obtained by Bhat- 
tal et al. (1998) were corrupted, both by the very large grav- 
ity smoothing length, and by the fact that the Jeans con- 
dition was violated (i.e. the Jeans mass was not resolved at 
all times). The inferences which they made are therefore not 
reliable. 



6 CONCLUSIONS 

We have conducted a series of simulations of star formation 
triggered by low- velocity collisions between low-mass molec- 
ular clumps. By implementing On-the-Fly Particle Splitting 
in these simulations, we have ensured that the Jeans con- 
dition is satisfied throughout the computational domain, at 
all times. Consequently we are confident that our simula- 
tions have not been corrupted by artificial fragmentation, 
and that real fragmentation has not been inhibited by in- 
adequate resolution. We are therefore now in a position to 
apply our SPH code with Particle Splitting to simulations 
of clump-clump collisions with more realistic initial condi- 
tions, for example clumps of unequal mass, clumps with in- 
ternal turbulence, and/or rotating clumps. The use of Par- 
ticle Splitting will afford significant advantages in terms of 
computational efficiency. 

In order to confirm the usefulness and reliability of this 
technique, we have repeated the Mo = 10M©, 6 = 0.2, 
M. — 10 simulation, using a standard code without Par- 
ticle Splitting, but introducing less massive (13 times less 
massive) particles, from the outset, and 13 times as many of 
them (i.e. 390,000 particles in total). The results differ only 
in small details, not in the overall statistics of fragmenta- 
tion, such as the number and mass of the protostars formed. 
Moreover, the standard simulation with 390,000 particles 
uses ~ 6 times more memory, and ~ 2 times more CPU 
than the simulation using On-the-Fly Particle Splitting with 
only 30,000 particles initially. 
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The simulations presented here demonstrate that the 
primary mode creating protostars, following a clump-clump 
collision, entails the formation of a shock compressed layer 
(while the gas is cooling radiatively) , its gravitational frag- 
mentation into one or more filaments (the number depend- 
ing on the lateral extent of the layer and the fragmentation 
scale), the break up of the filaments into cores, the conden- 
sation of primary protostars out of the cores, and the sub- 
sequent growth of the protostars by accretion along the fil- 
aments. Similar filamentary structures and cores have been 
found in the simulations of Pongracic et al. (1992), Klessen 
& Burkcrt (2000; 2001), Klessen (2001), Bate et al. (2002a; 
2002b; 2003), Bate & Bonnell (2005), Jappsen et al. (2005), 
Martel, Evans & Shapiro (2006). 

Increasing the clump mass Mo (with b and M held con- 
stant) increases the number of primary protostars, because 
the mass and lateral extent of the shock compressed layer 
is greater, and the fragmentation scale Lf rag is only moder- 
ately larger. Increasing b (with Mq and M held constant) 
reduces the number of primary protostars, because the lat- 
eral extent of the shock compressed layer is reduced (and 
the fragmentation length is increased somewhat by shear). 
Increasing b also increases the orbital angular momentum of 
the resulting ensemble of protostars. Increasing M (with Mo 
and b held constant) increases the number of primary pro- 
tostars because it reduces Lf rag . It also increases the orbital 
angular momentum of the resulting ensemble of protostars. 

For collisions with b > 0.6, no significant shock com- 
pressed layer forms, and the clump-clump collision does not 
trigger star formation. For collisions with b < 0.5, the fila- 
ments and cores that form have densities «h 2 > 10° cm" . 
Therefore it should be possible to map them in NH3 or 
CS, in nearby star formation regions. The protostars which 
condense out of a core have much higher central densities 
nu 2 > 10 10 cm~ 3 . Material accretes onto the protostars at 
rates in the range 1 to 5 x 10 _j Mq yr _1 , and this lasts for 
a time interval in the range 1 to 3 x 10 4 yrs. Therefore we 
should identify these protostars with the Class phase of 
evolution (Lada 1999). 

Hartmann (2002) has recently pointed out that young 
stars in Taurus form primarily within gaseous filaments, via 
fragmentation of the filaments; that the mean separation be- 
tween the young stars in Taurus is of order the local Jeans 
length; and that the protostellar cores are elongated along 
the filaments. This is very similar to the phenomenology 
which we observe in our simulations. In our simulations, the 
shock compressed layer created by the clump-clump collision 
fragments into filaments, and the filaments then fragment 
into cores; the separations between cores are of order the 
local Jeans length; and the individual cores are prolate and 
elongated parallel to the filaments. Our findings on filament 
fragmentation derived from clump-clump collision simula- 
tions agree with those of Jappsen et al. (2005) obtained from 
simulations of the turbulent ISM. They also fit the theoret- 
ical predictions of Larson (2005) on the fragmentation of 
filaments and the influence of thermodynamics. 

The prolateness of the cores in our simulations is only 
evident from the grey-scale plots we have presented, if one 
takes into account that molecular-line observations tend to 
be most sensitive to gas at the corresponding critical den- 
sity (~ 10 3 cm~ 3 to ~ 10 5 cm~ 3 , for the lines used to map 
cores in Taurus), and submillimeter observations are opti- 



cally thick for column-densities exceeding 4 x 10 25 H2 cm -2 . 
In other words one must look at the grey contours to see the 
structures which would be revealed by molecular-line or sub- 
millimeter continuum observations. Therefore, whereas the 
protostellar discs formed in our simulations present a range 
of shapes (dependent on whether they are viewed close to 
edge-on or close to pole-on), the cores within which these 
discs are embedded are almost invariably prolate and elon- 
gated along the filament from which they are condensing. 
The elongation along the filament arises because the cores 
are usually being fed by two opposing accretion streams from 
along the filament. 

Because the filaments are tumbling, the accretion flow 
onto a protostar tends to deliver ever increasing specific an- 
gular momentum, and therefore the protostar is spun up. 
This may lead to rotational instabilities which create sec- 
ondary companions to the primary protostars. However, the 
simulations have been terminated (due to limited compu- 
tational resource) before there is any clear evidence of this 
happening. 

Additionally, protostars condensing out of the same fil- 
ament will often fall together, and therefore they may inter- 
act. Ultimately, this could result in a merger, or in capture 
to form a binary, or in further fragmentation to produce 
secondary protostars, but again we have not followed the 
simulations for that long. The temptation is to introduce 
sink particles, but then the interaction is not properly mod- 
elled, because of the over-simplifications associated with the 
introduction of and the interaction between sink particles. 
Nevertheless, we have repeated one of the low-mass clump 
collision simulations using sink particles. This attempt has 
shown that, indeed, pairs of primary protostars can form 
close binaries through capture, as they move towards each 
other along a tumbling filament. It has also been shown that 
such binaries can be attended by massive circumbinary discs 
and that the formation of secondary companions to such a 
binary is possible through the interaction of its disc with the 
accretion flows onto it from along the filament. 

The SFE of low-mass clump collisions is estimated to be 
of order 20-30% based on the result of the simulation with 
sinks. We have asserted that the SFE decreases with in- 
creasing mass and, in particular, that the SFE of the 75M Q 
collisions is of order 10-15%. Both these values are in accor- 
dance with the findings of Hunter et al. (1986) from simu- 
lations of colliding gas flows. The SFEs estimated from our 
simulations are also consistent with those (20-30%) observed 
in a number of molecular clouds (Rengarajan 1984), which 
implies that clump-clump collisions can account, at least 
partly, for the star formation in those molecular clouds that 
can provide this triggering mechanism. 

Taken at face value, the above (theoretical and obser- 
vational) estimates for the SFE could be interpreted in the 
following sense: in the absence of any other triggering mech- 
anism, star formation in the ISM can be attributed mainly 
to low-mass clump collisions and not to collisions between 
higher-mass clumps. Evidently, this is not true in nature: 
first, these SFE estimates are only meaningful in a statistical 
sense, i.e. there is great variation from region to region; sec- 
ond, there is a number of complementary and/or mutually 
excluding mechanisms that are simultaneously in play. Col- 
lisions between low-mass clumps (and/or small-scale shocks 
produced in the turbulent ISM as part of larger-scale collid- 
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ing flows) appear to be one such mechanism able to explain 
part of the star formation observed in the Galaxy. 

Further progress will probably also necessitate simula- 
tions at even higher resolution, e.g. by employing multiple 
levels of particle splitting and using smaller sink particles 
(Kitsionas et al, in prep.). The detailed study of discs will 
also require the introduction of a more sophisticated algo- 
rithm to regulate the shear viscosity, for instance by us- 
ing the time-dependent formulation of Morris & Monaghan 
(1997) and/or the Balsara (1995) switch, the implementa- 
tion of a realistic continuity equation advocated by Imaeda 
& Inutsuka (2002) (despite the recent criticism by Monaghan 
(2006) of the need of the latter formalism) , and possibly all 
these should be attempted in combination with a different 
method for the velocity calculation, such as XSPH (Mon- 
aghan 2002). 
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